
# Optimal number of clusters according to silhouette widths
# (Rousseeuw quality index)
# *********************************************************

# Plot average silhouette widths (using Ward clustering) for all partitions 
# except for the trivial partition in a single group (k=1)
# First, create an empty vector in which the asw values will be written

grpdist <- function(X)
  {
  require(cluster)
  gr <- as.data.frame(as.factor(X))
  distgr <- daisy(gr, "gower")
  distgr
  }

silhouetteWidth <- function(x, clust, dissimilarity, type)
{
	if (type == 'silhouette')
	{
		asw <- numeric(nrow(x))
		for (i in 2:(nrow(x)-1))
		{
			sil <- silhouette(cutree(clust, k=i), dissimilarity)
			asw[i] <- summary(sil)$avg.width
		}
		k.best <- which.max(asw)
		print(k.best)
		plot(1:nrow(x), asw, type="h", 
			 main="Silhouette-optimal number of clusters", 
			 xlab="k (number of groups)", ylab="Average silhouette width")
			 axis(1, k.best, paste("optimum",k.best,sep="\n"), col="red", font=2,
			 col.axis="red")
		points(k.best, max(asw), pch=16, col="red", cex=1.5)
	}
	else if (type == 'mantel')
	{
		kt <- data.frame(k=1:nrow(x), r=0)
		for (i in 2:(nrow(x)-1))
		{
			gr <- cutree(clust, i)
			distgr <- grpdist(gr)
			mt <- cor(dissimilarity, distgr, method="pearson")
			kt[i,2] <- mt
		}
		k.best <- which.max(kt$r)
		plot(kt$k, kt$r, type="h", main="Mantel-optimal number of clusters", 
			 xlab="k (number of groups)", ylab="Pearson's correlation")
			 axis(1, k.best, paste("optimum", k.best, sep="\n"), col="red", font=2,
			 col.axis="red")
		points(k.best, max(kt$r), pch=16, col="red", cex=1.5)
	}
}

# Silhouette plot of the final partition
# **************************************
silhouetteFinal <- function(x, clust, dissimilarity, numClust)
{
	cutg <- cutree(clust, k=numClust)
	sil <- silhouette(cutg, dissimilarity)
	silo <- sortSilhouette(sil)
	rownames(silo) <- row.names(x)[attr(silo,"iOrd")]
	par(mar=c(8,5,5,5))
	plot(silo, main="Silhouette plot", cex.names=0.8, col=cutg+1, nmax.lab=100)
}
